import string
import copy
import scipy
import Tkinter, tkFileDialog
import numpy as np
import pandas as pd
from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt
import os
import sys
import cPickle
import matplotlib.units
sys.path.append(os.path.abspath("C:\Users\Scherer Lab E\Documents\GitHub\Python_Data_Analysis"))
from common_functions import *
%matplotlib inline
#%config InlineBackend.figure_format = 'retina'
#import seaborn as sns

os.chdir("C:\Users\Scherer Lab E\Box Sync\Ring\Figure Data and Scripts\Data")

#matplotlib.rcParams

matplotlib.rcParams['font.serif']

[u'Bitstream Vera Serif',
 u'DejaVu Serif',
 u'New Century Schoolbook',
 u'Century Schoolbook L',
 u'Utopia',
 u'ITC Bookman',
 u'Bookman',
 u'Nimbus Roman No9 L',
 u'Times New Roman',
 u'Times',
 u'Palatino',
 u'Charter',
 u'serif']

#matplotlib.rcParams['font.family'] = 'Vera'
matplotlib.rcParams['figure.dpi'] = 200.0
matplotlib.rcParams['axes.labelsize'] = 15
matplotlib.rcParams['axes.titlesize'] = 18
matplotlib.rcParams['xtick.labelsize'] = 'large'
matplotlib.rcParams['ytick.labelsize'] = 'large'
matplotlib.rcParams['mathtext.fontset'] = 'stix'
matplotlib.rcParams['font.family'] = 'STIXGeneral'
matplotlib.rcParams['legend.handletextpad'] = 0.1
matplotlib.rcParams['legend.numpoints'] = 1
matplotlib.rcParams['legend.scatterpoints'] = 1
matplotlib.rcParams['legend.borderaxespad'] = 0.9
#matplotlib.rcParams['mathtext.default'] = 'regular'
#matplotlib.rcParams['font.family'] = 'serif'
#matplotlib.rcParams['font.serif'] = 'Times'

#matplotlib.rcParams['text.usetex'] = True
#matplotlib.rc('text.latex', preamble=r'\usepackage{cmbright}')

def add_axes_label_inches(ax, (right, down), string, **kwargs):
    """A helper function to add a text lable a specific distance from the top left corner in inches
    
    This function was created so that when constructing multi-panel figures for a manuscript all the axes labels
    (e.g. a, b, c, d...) are a consistent distance from the corner of the axis. The relative coordinates used for text
    labels can result in a label being a different physical distance compared to subplots of different sizes.
    
    :param fig: A mpl.figure object that contains the axis you want to add text to
    :param ax: A mpl.axes object that you want to add text to
    :param (right,down): Distance in inches you want the text label to be from the top left corner of the axis
    :param string: The text string you want displayed at that location."""
    fig = ax.get_figure()
    fig_size = fig.get_size_inches()
    ax_bbox = ax.get_position()
    ax_rect_inches = ax_bbox.x0*fig_size[0], ax_bbox.y0*fig_size[1], ax_bbox.x1*fig_size[0], ax_bbox.y1*fig_size[1]
    text_location_inches = [right, ax_rect_inches[3]-ax_rect_inches[1]-down]
    text_position_rel_coors = text_location_inches[0]/(ax_rect_inches[2]-ax_rect_inches[0]), text_location_inches[1]/(ax_rect_inches[3]-ax_rect_inches[1])
    return ax.text(text_position_rel_coors[0], text_position_rel_coors[1], string, transform=ax.transAxes, va='top', ha='left', **kwargs)

fig = plt.figure(figsize=[3,4])
ax = plt.subplot()
add_axes_label_inches(ax, (0.3,0.3), 'a', fontsize=14)
plt.close()

Figure 5

Velocity of Particles around Ring with Force Maps

Exerimental Velocity around Ring

Load data from Mov_11121401 the single particle ramping over the nanoplate (see Ana_16011101). Only the L=4 part of the trajectory is used because it looks the best

import cPickle
f = open("C:\Users\Scherer Lab E\Box Sync\Ring\Figure Data and Scripts\Data\Fig_7.pkl")
trajs_pl = cPickle.load(f)
f.close()

def avg_st_dev_velocity_vs_theta(data_frame, t, mag_slider=1.6, theta_bin_num=18, ax=None):
    um_conv=6.5/60/mag_slider/2
    copy = data_frame.drop_duplicates(subset=['frame','track id']).copy()
    
    
    dist = lambda x: (((x.shift(-t)['x pos']-x.shift(t)['x pos'])**2 + (x.shift(-t)['y pos']-x.shift(t)['y pos'])**2)**(0.5))/2.0
    grouped = copy.groupby(['track id'])
    
    # Need to stack if only one group
    if len(grouped) == 1:
        result = copy.groupby(['track id']).apply(dist)
        result = result.stack()
    else:    
        result = copy.groupby(['track id']).apply(dist)
    
    copy['displacement'] = result.reset_index(level='track id').drop('track id', axis=1)
    copy = copy.dropna(axis=0, subset=['displacement']).copy()
    
    copy['velocity'] = copy['displacement']*um_conv*90.0
    
    theta_bins = np.linspace(0, 360, num=theta_bin_num)
    theta_bin_mid_points = (np.roll(theta_bins,-1)[:-1] + theta_bins[:-1])/2.0
    copy['theta_bin'] = pd.cut(copy['theta'], theta_bins)
    avgs = copy[['velocity','theta_bin']].groupby('theta_bin').mean().values[:,0]
    std_dev = copy[['velocity','theta_bin']].groupby('theta_bin').std().values[:,0]
    avgs = np.roll(avgs, -5)
    std_dev = np.roll(std_dev, -5)    
    ax.plot(theta_bin_mid_points*np.pi/180.0, avgs, '-', color="#4C72B0", mec="#4C72B0")
    ax.fill_between(theta_bin_mid_points*np.pi/180.0, avgs+std_dev, avgs-std_dev, facecolor="#4C72B0", alpha=0.5, lw=0.0)
    #plt.plot(theta_bin_mid_points, avgs-std_dev)
    plt.ylabel('Velocity ($\mathregular{\mu}$m/sec)')
    x_ticks = np.array([0,0.5*np.pi,1*np.pi,1.5*np.pi,2*np.pi])
    x_label = [r"$0$", r"$\pi/2$", r"$\pi$", r"$3\pi/2$", r"$2\pi$"]
    plt.xticks(x_ticks, x_label)
    #plt.xlabel('Theta')

%matplotlib inline
ax = plt.subplot(projection='polar')
avg_st_dev_velocity_vs_theta(trajs_pl[9], 1, theta_bin_num=20, ax=ax)
#plt.close()
plt.show()

Simulation of Velocity around Ring

velocity_position_L4.mat is the velocity of a single particle in an L=4 trap. From Nishant:

This file contains the velocity of a single particle in an L=4 trap. theta_l4 (radians) and vmag_l4 (m/s) are the raw velocity and position data. mean_vmag is the binned average and theta_bin is the bin location. std_dev is the standard deviation for the binned velocity data.

theta-bin is actually the bin mid points for the mean velocity values (mean_vmag)

'''Load the simulation .mat file'''

import scipy.io as sio
vel_sim_dat = sio.loadmat('velocity_position_L4.mat')

# A list of useless keys to skip in each .mat file
skip_list = ['__version__', '__header__', '__globals__']

# Turn all the arrays into 1D arrays for the ang_msd_fit_dat
new_vel_sim_dat = {}
for key, val in vel_sim_dat.iteritems():
    if key in skip_list:
        continue
    #print val.shape, key
    if key in ['theta_bin']:
        new_vel_sim_dat[key] = val[0, :]
    else:
        new_vel_sim_dat[key] = val[:, 0]
vel_sim_dat = new_vel_sim_dat

fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})
theta = vel_sim_dat['theta_bin']
mean_vel = vel_sim_dat['mean_vmag']*10**6
std_dev = vel_sim_dat['std_dev']*10**6
ax.plot(theta, mean_vel, '-o', color="#4C72B0", mec="#4C72B0")
ax.fill_between(theta, mean_vel+std_dev, mean_vel-std_dev, facecolor="#4C72B0", alpha=0.5, lw=0.0)
x_ticks = np.array([0,0.5*np.pi,1*np.pi,1.5*np.pi,2*np.pi])
x_label = [r"$0$", r"$\frac{\pi}{2}$", r"$\pi$", r"$3\pi/2$", r"$2\pi$"]
plt.xticks(x_ticks, x_label)
plt.ylabel('Velocity ($\mathregular{\mu}$m/sec)')
plt.xlabel('Theta')
plt.title("Simulations Over Plate L=4")

#plt.close()
plt.show()

Import Simulation Force Map Data

The force map and electric field intensity data are in forcemap_l4.mat are for L=4 from the simulations. From Nishant:

This file contains the optical force vectors and positions, and the electric field intensity in the ring trap. force_x (force_y) is the x (y) component of the force. pos_x (pos_y) is the position along x (y) axis. field_intensity is the electric field intensity and x, y are just x and y axes for the field intensity. The idea here is to plot the field intensity as the background and make a quiver plot with the force vectors.

'''Load the simulation .mat file'''

import scipy.io as sio
force_sim_dat = sio.loadmat('forcemap_l4.mat')

# A list of useless keys to skip in each .mat file
skip_list = ['__version__', '__header__', '__globals__']

# Turn all the arrays into 1D arrays for the ang_msd_fit_dat
new_force_sim_dat = {}
for key, val in force_sim_dat.iteritems():
    if key in skip_list:
        continue
    #print val.shape, key
    if key in ['x', 'y']:
        new_force_sim_dat[key] = val[0, :]
    else:
        new_force_sim_dat[key] = val
force_sim_dat = new_force_sim_dat

def add_two_polar_coords(ax, aspect='equal'):
    """A helper function that adds 0pi and pi/2 axis labels on a twinx and
    twin y axis.
    
    :param ax: A mpl.axes object you want to add polar labels to
    :param aspect: The aspect ratio you want the final plot to have
    """
    ax_y_polar = ax.twinx()
    ax_y_polar.set_ylim(ax.get_ylim())
    ax_y_polar.set_yticks([0])
    ax_y_polar.set_yticklabels(['$0\pi$'])
    
    ax_x_polar = ax.twiny()
    ax_x_polar.set_xlim(ax.get_xlim())
    ax_x_polar.set_xticks([0])
    ax_x_polar.set_xticklabels(['$\pi/2$'])
    
    ax.set(adjustable='box-forced', aspect=aspect)
    ax_y_polar.set(adjustable='box-forced', aspect=aspect)
    ax_x_polar.set(adjustable='box-forced', aspect=aspect)

def add_four_polar_coords(ax):
    """A helper function that adds 0pi, pi/2, pi, and 3pi/2 axis labels 
    on the inside of an axis following the unit circle.
    
    :param ax: A mpl.axes object you want to add polar labels to
    """
    # Add the polar labels
    ax.text(0.97, 0.5, '$\mathbf{0}$', fontdict=dict(color='white', size='large'),
         weight='extra bold', transform=ax.transAxes, va='center', ha='right')
    ax.text(0.5, 0.97, '$\mathbf{\pi/2}$', fontdict=dict(color='white', size='large'),
             weight='extra bold', transform=ax.transAxes, va='top', ha='center')
    ax.text(0.03, 0.5, '$\mathbf{\pi}$', fontdict=dict(color='white', size='large'),
             weight='extra bold', transform=ax.transAxes, va='center', ha='left')
    ax.text(0.5, 0.03, '$\mathbf{3\pi/2}$', fontdict=dict(color='white', size='large'),
             weight='extra bold', transform=ax.transAxes, va='bottom', ha='center')
    
    # Change ticks in graph
    xticklines = ax.get_xticklines()
    yticklines = ax.get_yticklines()
    
    xcenter = len(xticklines)/2
    ycenter = len(yticklines)/2
    
    xticklines[xcenter].set(color='white', markeredgewidth=2, markersize=5, zorder=0)
    xticklines[xcenter-1].set(color='white', markeredgewidth=2, markersize=5)
    yticklines[ycenter].set(color='white', markeredgewidth=2, markersize=5)
    yticklines[ycenter-1].set(color='white', markeredgewidth=2, markersize=5)

%matplotlib inline

force_fig = plt.figure(figsize=(4,4))
# Plot electric field intensity
field_x = force_sim_dat['x']/10**3
field_y = force_sim_dat['y']/10**3
field_int = force_sim_dat['field_intensity']
plt.pcolormesh(field_x, field_y, field_int, cmap='viridis')
plt.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])

# Plot the force vectors
force_pos_x = force_sim_dat['pos_x']/10**3
force_pos_y = force_sim_dat['pos_y']/10**3
force_mag_x = force_sim_dat['force_x']/10**3
force_mag_y = force_sim_dat['force_y']/10**3
plt.quiver(force_pos_x[::2, ::2], force_pos_y[::2, ::2], force_mag_x[::2, ::2], force_mag_y[::2, ::2], color="#4D0E10")
#plt.quiver(force_pos_x[:, :], force_pos_y[:, :], force_mag_x[:, :], force_mag_y[:, :], color="#4D0E10")

plt.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])

add_two_polar_coords(plt.gca())

plt.close()

Import the Force Map data from the Simulation L=0

This data is from Nishant's simulation for L=0. From Nishant:

File: forcemap_l0.mat x, y, field_intensity variables are the same as above. pos_x, pos_y and force_x, force_y are the coordinates and components of the force vectors for L=0.

'''Load the simulation .mat file'''

import scipy.io as sio
force_sim_dat_L0 = sio.loadmat('forcemap_l0.mat')

# A list of useless keys to skip in each .mat file
skip_list = ['__version__', '__header__', '__globals__']

# Turn all the arrays into 1D arrays for the ang_msd_fit_dat
new_force_sim_dat_L0 = {}
for key, val in force_sim_dat_L0.iteritems():
    if key in skip_list:
        continue
    #print val.shape, key
    if key in ['x', 'y']:
        new_force_sim_dat_L0[key] = val[0, :]
    else:
        new_force_sim_dat_L0[key] = val
force_sim_dat_L0 = new_force_sim_dat_L0

# Plot the Force Map
force_fig = plt.figure(figsize=(4,4))
# Plot electric field intensity
field_x = force_sim_dat_L0['x']*10**6
field_y = force_sim_dat_L0['y']*10**6
field_int = force_sim_dat_L0['field_intensity']
plt.pcolormesh(field_x, field_y, field_int, cmap='viridis')
plt.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])

# Plot the force vectors
force_pos_x = force_sim_dat_L0['pos_x']*10**6
force_pos_y = force_sim_dat_L0['pos_y']*10**6
force_mag_x = force_sim_dat_L0['force_x']*10**6
force_mag_y = force_sim_dat_L0['force_y']*10**6
plt.quiver(force_pos_x[::2, ::2], force_pos_y[::2, ::2], force_mag_x[::2, ::2], force_mag_y[::2, ::2], color="#4D0E10")

plt.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])

add_two_polar_coords(plt.gca())

plt.close()

Combining the Plots

from matplotlib.gridspec import GridSpec

vel_exp_and_sim_force = plt.figure(figsize=(12,4))

gs = GridSpec(2,3, width_ratios=[1,1.5,1.5])


ax1 = plt.subplot(gs[0,0])
avg_st_dev_velocity_vs_theta(trajs_pl[9], 1, theta_bin_num=20, ax=ax1)
#plt.title("Experiment L=4")
plt.setp(ax1.get_xticklabels(), visible=False)
ax1.set_ylim(ymax=1.2*ax1.get_ylim()[1])
ax1.text(0.5, .94, 'Experiment L=4', transform=ax1.transAxes,
         va='top', ha='center', fontsize=12)

ax2 = plt.subplot(gs[1,0])
theta = vel_sim_dat['theta_bin']
mean_vel = vel_sim_dat['mean_vmag']*10**6
std_dev = vel_sim_dat['std_dev']*10**6
ax2.plot(theta/np.pi, mean_vel, '-', color="#4C72B0", mec="#4C72B0")
ax2.fill_between(theta/np.pi, mean_vel+std_dev, mean_vel-std_dev, facecolor="#4C72B0", alpha=0.5, lw=0.0)
x_ticks = np.array([0,0.5,1,1.5,2])
x_label = [r"$0$", r"$\pi/2$", r"$\pi$", r"$3\pi/2$", r"$2\pi$"]
plt.xticks(x_ticks, x_label)
plt.ylabel('Velocity ($\mathregular{\mu}$m/sec)')
plt.xlabel(r'$\theta$')
#plt.title("Simulation L=4")
ax2.set_ylim(ymax=1.2*ax2.get_ylim()[1])
ax2.text(0.5, .94, 'Simulation L=4', transform=ax2.transAxes,
         va='top', ha='center', fontsize=12)

# Plot the Force Map L=0
ax3 = plt.subplot(gs[:,1])
# Plot electric field intensity
field_x = force_sim_dat_L0['x']*10**6
field_y = force_sim_dat_L0['y']*10**6
field_int = force_sim_dat_L0['field_intensity']
plt.pcolormesh(field_x, field_y, field_int, cmap='viridis')
plt.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])

# Plot the force vectors
force_pos_x = force_sim_dat_L0['pos_x']*10**6
force_pos_y = force_sim_dat_L0['pos_y']*10**6
force_mag_x = force_sim_dat_L0['force_x']*10**6
force_mag_y = force_sim_dat_L0['force_y']*10**6
plt.quiver(force_pos_x[::2, ::2], force_pos_y[::2, ::2], force_mag_x[::2, ::2], force_mag_y[::2, ::2], color="#4D0E10")
plt.xlabel('x ($\mathregular{\mu}$m)')
plt.ylabel('y ($\mathregular{\mu}$m)')
plt.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])

#add_two_polar_coords(plt.gca())
ax3.set(adjustable='box-forced', aspect='equal')
add_four_polar_coords(ax3)

# Plot Force Map L=4
ax4 = plt.subplot(gs[:,2])
# Plot electric field intensity
field_x = force_sim_dat['x']/10**3
field_y = force_sim_dat['y']/10**3
field_int = force_sim_dat['field_intensity']
plt.pcolormesh(field_x, field_y, field_int, cmap='viridis')
plt.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])

# Plot the force vectors
force_pos_x = force_sim_dat['pos_x']/10**3
force_pos_y = force_sim_dat['pos_y']/10**3
force_mag_x = force_sim_dat['force_x']/10**3
force_mag_y = force_sim_dat['force_y']/10**3
ax4.quiver(force_pos_x[::2, ::2], force_pos_y[::2, ::2], force_mag_x[::2, ::2], force_mag_y[::2, ::2], color="#4D0E10")
plt.axis([field_x.min(), field_x.max(), field_y.min(), field_y.max()])

plt.xlabel('x ($\mathregular{\mu}$m)')
plt.ylabel('y ($\mathregular{\mu}$m)')

#add_two_polar_coords(ax4)
ax4.set(adjustable='box-forced', aspect='equal')
add_four_polar_coords(ax4)

plt.tight_layout()

gs.update(hspace=0.1)

add_axes_label_inches(ax1, (0.075,0.075), 'a', fontsize=14, weight='extra bold')
add_axes_label_inches(ax2, (0.075,0.075), 'b', fontsize=14, weight='extra bold')
text_c = add_axes_label_inches(ax3, (0.075,0.075), 'c', fontsize=14, weight='extra bold', color='white')
text_d = add_axes_label_inches(ax4, (0.075,0.075), 'd', fontsize=14, weight='extra bold', color='white')
ax3.annotate('', xy=(1,17), xycoords='axes points', 
            xytext=(24,0), textcoords='offset points', 
            arrowprops=dict(edgecolor='red', facecolor='red', arrowstyle='<|-|>',
                            connectionstyle='arc3,rad=0'))
ax3.annotate('E', xy=(0,0), xytext=(13,11), textcoords='axes points',
            horizontalalignment='center', verticalalignment='center', color='red', weight='extra bold')
ax4.annotate('', xy=(1,17), xycoords='axes points', 
            xytext=(24,0), textcoords='offset points', 
            arrowprops=dict(edgecolor='red', facecolor='red', arrowstyle='<|-|>',
                            connectionstyle='arc3,rad=0'))
ax4.annotate('E', xy=(0,0), xytext=(13,11), textcoords='axes points',
            horizontalalignment='center', verticalalignment='center', color='red', weight='extra bold')

<matplotlib.text.Annotation at 0x285add68>

save_dir = "C:\Users\Scherer Lab E\Box Sync\Ring\Figure Data and Scripts\Figures"
vel_exp_and_sim_force.savefig(save_dir+"\Fig_5.png", dpi=800)